function [g,h]=SQPcon(x)
% output table like shown in class
% convert to constrain function

% x design variable vector, 2 stage compound gear train
% x(1) N1 - # of teeth input driving gear (stage 1)
% x(2) N2 - # of teeth input driven gear (stage 1)
% x(3) N3 - # of teeth output driving gear (stage 2)
% x(4) N4 - # of teeth output driven gear (stage 2)
% x(5) P12 - diametral pitch of gears 1/2
% x(6) P34 - diametral pitch of  gears 3/4

%%% VARIABLES
phi = 20*pi/180; % pressure angle [rad], increase to reduce undercutting

switch '2' % 1 for only first 4 constraints to find a design space, 2 full analysis
    case '1'
        g(:,1) = (x(:,1)+x(:,2))./x(:,5)==(x(:,3)+x(:,4))./x(:,6); % check centerline constraint
        g(:,2) = (x(:,1)./x(:,2)).*(x(:,3)./x(:,4)) >= 17; % minimum overall gear ratio
        g(:,3) = (x(:,1)./x(:,2)).*(x(:,3)./x(:,4)) <= 21; % maximum overall gear ratio
        return
    case '2'
        % just keep goin
end

N1 = x(:,1); 
N2 = x(:,2);
N3 = x(:,3);
N4 = x(:,4);
P12 = x(:,5);
P34 = x(:,6);

% split variables with a column for each gear to make some calculations
% easier

% to access individual variables using this format:
%gear 1 VAR(:,1)
%pinion 2 VAR(:,2)
%gear 3 VAR(:,3)
%pinion 4 VAR(:,4)

P(:,1) = P12;
P(:,2) = P12;
P(:,3) = P34;
P(:,4) = P34;
N(:,1) = N1;
N(:,2) = N2;
N(:,3) = N3;
N(:,4) = N4;

GR = (N1./N2).*(N3./N4); %overall gear ratio

%%% UNDERCUT ANALYSIS 
%(pg 686)- used for constraints eventually
k = 1.0; % full teeth
% k = 0.8; % stub teeth
Mg = (N1./N2); % for input gear set
Np1 = ceil((2*k)*(Mg+sqrt(Mg.^2+(1+2*Mg)*sin(phi)^2))./((1+2*Mg)*sin(phi)^2));
Mg = (N3./N4); % for output gear set
Np2 = ceil((2*k)*(Mg+sqrt(Mg.^2+(1+2*Mg)*sin(phi)^2))./((1+2*Mg)*sin(phi)^2));



%%% CONTACT RATIO
d = N./P; % pitch diameter [inch] (eq 13-1, pg 676)
r = d/2; % pitch radius [inch]
pc(:,1) = pi./P(:,1); % circular pitch, set 1 (eq 13-4, pg 676)
pc(:,2) = pi./P(:,3); % circular pitch, set 2
rb = r*cos(phi); % base circle radius [inch] (eq 13-6, pg 680)
a = 1./P; % addendum, (pg 680)

LAB(:,1) = sqrt((r(:,2)+a(:,1)).^2-rb(:,2).^2) + sqrt((r(:,1)+a(:,1)).^2-rb(:,1).^2) - (r(:,2)+r(:,1))*sin(phi); % line of action, set 1, (eq 14-25, pg 755)
%LAB(:,2) = LAB(:,1);
LAB(:,2) = sqrt((r(:,4)+a(:,3)).^2-rb(:,4).^2) + sqrt((r(:,3)+a(:,3)).^2-rb(:,3).^2) - (r(:,4)+r(:,3))*sin(phi); % line of action, set 2
%LAB(:,4) = LAB(:,3);  

mct = LAB./(pc*cos(phi)); % contact ratio, (eq 13-9, pg 685)
mc = min(mct,[],2); % minimum contact ratio

%%% FORCE ANALYSIS
T = 38.6; % torque [ft*lb] from problem statement
% Wt = shaft torque/pitch circle radius
Wt(:,1) = T./(r(:,1)/12); % tangential load on gear 1 [lb]
Wt(:,2) = Wt(:,1); % tangential load on pinion 2 [lb], contact force between gear set must be equal
Wt(:,3) = Wt(:,2).*(r(:,2)/12)./(r(:,3)/12); % tangential load on gear 3 [lb]
Wt(:,4) = Wt(:,3); % tangential load on pinion 4 [lb], contact force between gear set must be equal

% Face Width
% initial approximation ROT pg. 776
F = 4*pi./P; % [in]

% AGMA Stress Correction Factors
Ko = 1.25; % Overload factor p.758, figs(p.766), Starkey: 1-1.25

V12 = 200*(pi*d(:,1)/12); %pitchline velocity gear set 1 [ft/min] (eq. 13-34, pg 707)
V34 = 200*(N1./N2).*(pi*d(:,3)/12); %pitchline velocity gear set 2 [ft/min] (eq. 13-34, pg 707)
Qv = 8; % Starkey Recommends Qv of 7,8,9, will pick 8 initially
B = 0.25*(12-Qv)^(2/3); % kv calculation constant p756
A = 50+56*(1-B); % kv calculation constant p756
Vt_max = (A+(Qv-3))^2; % max pitchline velocty [ft/min] pg 756
Kv12 = ((A+sqrt(V12))/A).^B; % dynamic factor, V in [ft/min]
Kv34 = ((A+sqrt(V34))/A).^B; % dynamic factor, V in [ft/min]
Kv= [Kv12,Kv12,Kv34,Kv34];

Y = zeros(size(Kv));
J = Y;
J(:,1) = 0.5; % lewis geometry factor p.753 (J in book), ~ 125 teeth
J(:,2) = 0.35; % lewis geometry factor p.753 (J in book) ~20 teeth
J(:,3) = 0.5; % lewis geometry factor p.753 (J in book) ~125 teeth
J(:,4) = 0.35; % lewis geometry factor p.753 (J in book) ~20 teeth

m_n = 1; % pg 755, also used for I
Y(:,1) = 0.45; % lewis form factor (pg 738 table 14-2) for 100 teeth [ ]
Y(:,2) = 0.35; % for 28 teeth (conservative)
Y(:,3) = 0.45; % for 100 teeth
Y(:,4) = 0.35; % for 28 teeth (conservative

Ks = max(1.192*(F.*sqrt(Y)./P).^0.0535,1); % size factor p.759

Cmc = 1; % uncrowned teeth pg. 760
Cpf = F./(10*d) - .0375 + .0125*F; %  Assume 1 < F < 17, pg 760
Cpm = zeros(size(Kv));
Cpm(:,1) = 1; % Gear 1 centered between bearings, pg 760
Cpm(:,2) = 1.1; % Gear 2 not centered btwn brngs, pg 760
Cpm(:,3) = 1.1; % Gear 3 not centered btwn brngs, pg 760
Cpm(:,4) = 1; % Gear 4 centered between bearings, pg 760
Cma = .127 + .0158*F - .93E-4*F.^2; % Commercial, enclosed units, pg 760
Ce = 1; % other conditions, pg 760
Km = 1 + Cmc*(Cpf.*Cpm+Cma.*Ce); % load distribution factor, eq 14-30 pg 759

Kb = 1; % rim-thickness factor (1 for full gears)

Cp = 2300; % steel on steel, pg 757, tbl 14-8 [psi^0.5]
Cf = 1; % no correlation established, just use 1 pg758

GG = ([N1./N2,N1./N2,N3./N4,N3./N4]); % define gear ratio for each gear (4 columns (1 for each gear) x (number of designs rows)
I = (cos(phi)*sin(phi)/(2*m_n)) * (GG./(GG+1)); %geometry factor, eq 14-23, pg 755

Sigma_Bend = Wt*Ko.*Kv.*Ks.*P.*Km*Kb./(F.*J); % AGMA Bending Stress, eq 14-15, pg 746, 766
Sigma_Cont = Cp*(Wt*Ko.*Kv.*Ks.*Km*Cf./(d.*F.*I)).^(1/2); %AGMA Contact Stress, eq 14-16, pg 746, 766

%AGMA Strength Correction Factors
%**NEED TO ADD HARDNESS INPUT**
Hb = 390; %%% arbitrary choice!
St = 77.3*Hb + 12800; % allowable bending stress, thru-hardened steel, assume grade 1, figure 14-2, pg 747 [psi]
Yn = 0.8; % Stress cycle factor, 10^10 cycles @ 3600 RPM = 5yr life.  Average of Ynmax = .864, Ynmin = .743. fig 14-14, pg 763
Kt = 1; % standard temperature p.764
Kr = 1; % Reliability factor, 99%, table 14-10, pg 764

Sc = 322*Hb + 29100; % allowable contact stress, thru-hardened steel, grade 1, fig 14-5, pg 750 [psi]
Zn = .845; % stress cycle life factor, 10^10 cycles @ 3600 RPM = 5 yr life.  Average of Znmax = .881, Znmin = .809. Fig 14-15, pg 763
Ch = 1; % hardness ratio factors, Assume hardness of gear and pinion are equal -> A' = 0.  eq 14-36, pg 761

% safety factors
Sf = St*Yn./(Kt*Kr*Sigma_Bend); %Bending Safety Factor, eq 14-41, pg 765
Sh = Sc*Zn*Ch./(Kt*Kr*Sigma_Cont); % Pitting Safety Factor, eq. 14-42, pg 765

Sigma_Bend_All = St*Yn./(Sf*Kt*Kr); %Bending Allowable Stress, eq. 14-17, pg 749
Sigma_Cont_All = Sc*Zn*Ch./(Sh*Kt*Kr); %Contact Allowable Stress Calculation, eq 14-18, pg 750 [kpsi]

%%% EVALUATE COSTRAINTS
g(:,1) = Np1/N2 -1; % check undercut gear set 1 
g(:,2) = Np2/N4 - 1; % check undercut gear set 2
g(:,3) = 17/GR -1; % minimum overall gear ratio
g(:,4) = GR/21 - 1; % maximum overall gear ratio
g(:,5) = 1.2/mc - 1; % contact ratio greater than 1.2
g(:,6:9) = 1.5./Sf - 1; % Bending Factor of Safety Gear INCOMPLETE
g(:,10:13) = 1.2./Sh - 1; % Contact Factor of Safety Gear INCOMPLETE

h = [(N1+N2)./P12 - (N3+N4)./P34]; % nonlineary equality constraints
